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Abstract 

The Gaussian random field (GRF) and the Gaussian Markov random field (GMRF) have 
been widely used to accommodate spatial dependence under the generalized linear mixed 
model framework. These models have limitations rooted in the symmetry and thin tail of the 
Gaussian distribution. We introduce a new class of random fields, termed transformed GRF 
(TGRF), and a new class of Markov random fields, termed transformed GMRF (TGMRF). 
They are constructed by transforming the margins of GRFs and GMRFs, respectively, to de- 
sired marginal distributions to accommodate asymmetry and heavy tail as needed in practice. 
The Gaussian copula that characterizes the dependence structure facilitates inferences and ap- 
plications in modeling spatial dependence. This construction leads to new models such as 
gamma or beta Markov fields with Gaussian copulas, which can be used to model Poisson 
intensity or Bernoulli rate in a spatial generalized linear mixed model. The method is natu- 
rally implemented in a Bayesian framework. We illustrate the utility of the methodology in an 
ecological application with spatial count data and spatial presence/absence data of some snail 
species, where the new models are shown to outperform the traditional spatial models. The 
validity of Bayesian inferences and model selection are assessed through simulation studies 
for both spatial Poisson regression and spatial Bernoulli regression. 

Some key words: Bayesian inference, beta field, gamma field, Gaussian copula, generalized linear mixed model 
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A Gaussian random field (GRF) is a stochastic process whose finite dimensional marginal distribu- 
tion of any dimension is Gaussian. A GRF that enjoys the Markov property is a Gaussian Markov 
random field (GMRF) and can be represented by an undirected graph. Specifically, a random vec- 
tor Z follows a GMRF with respect to some labeled graph Q if Z is a GRF with precision matrix 
Q such that ^ if and only if i and j are connected in graph Q ( |Rue and Held 2005 ). A 
GMRF model is natural when specification of precision Q is easier than specification of covari- 
ance, £ = Q 1 . From the Markov property, the conditional distribution of node i in Q given all the 
rest is fully specified by its neighbors. The precision matrix Q plays an essential role here and the 
sparsity of Q facilitates fast sampling algorithms, which is important in sampling based inferences 
(Rue| 2001[ ). GMRFs are analytically tractable, can be easily used in hierarchical models, and fit 
nicely in a Bayesian framework. 

The GMRF model has been widely used in a variety of fields. Along with the public concerns 
in the environment and public health, recent applications have surged in environmental sciences 
(e.g., |Wikle et a!4|1998||Huerta et aL||2004t|Rue et al4|2004| ) and epidemiology (e.g., |Besag et al.[ 
T99T| |Knorr-Held and Besag] [T998t |Held and Ruej |2002l |5dimid and Held} [20041 ). In particular, 



GMRFs are used as random effects to account for spatial dependence in the generalized linear 
mixed model (GLMM) framework; see Rue and Held ( |2005| ) and references therein. Popular as 
they are, GMRFs are limited in accommodating asymmetry or heavy tails in practice because of the 
marginal Gaussian property. To the best of our knowledge, Markov random fields with arbitrary 

21 margins are underdeveloped for modeling data with asymmetry or heavy tails. 

22 We propose random field that are transformed from a GRF and GMRF such that the marginal 

23 distributions are of any desired form, through the probability integral transformation and its in- 

24 verse. By Sklar's theorem (Sklar |1959[ ), any continuous multivariate distribution can be uniquely 

25 represented by its marginal distributions and a copula which characterizes the dependence struc- 
ture. Since copulas are invariant to monotonic transformations, the dependence structure of the 
new field, in terms of copula, remains the same as that of the GRF or GMRF. Such construction 
leads to new random fields with desired margins combined with Gaussian copulas, such as gamma 
fields or beta fields, which may be used to model spatial Poisson intensity or Bernoulli rate. The 
transformed GRF (TGRF) and the transformed GMRF (TGMRF) share some of the properties of 
the GRF and GMRF, respectively. Due to the marginal transformations, they allow a general and 
flexible representation that can easily accommodate asymmetry as well as heavy tail behavior that 
are often observed in empirical data. 

The rest of the article is organized as follows. In Section [2j TGRFs and TGMRFs are defined 
and some of their important properties are presented. In Section [3j TGMRFs are incorporated 
into a GLMM framework to accommodate spatial dependence. In Section [4] and Section [5} the 
proposed models are applied to the count data and presence/absence data, respectively, of some 
snail species in an ecological study. The statistical inferences are made in the Bayesian framework. 
The performance of the inferences and model selection in both applications are investigated in 
simulations that mimic the real data. A discussion concludes in Section |6j 
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1 2 A General Class of Random Fields 

2 For ease of notation, we present the definition of a TGRF and a TGMRF in the context of finite 

3 dimension n in the sequel. For random fields indexed by elements in some space, the definition 

4 applies to n-dimensional marginal distributions for any n. 

5 Suppose that e = (ei, . . . , e n ) T is n-dimensional standard multivariate normal with mean 

6 and correlation matrix, St, denoted as iV n (0, Define a random vector Z = (Z\, . . . , Z n ) T 

7 through 

Z i = Ff 1 {$(£,))}, i = l,...,n, (1) 

e where F { is the distribution function of an absolutely continuous variable and $ is the distribution 

9 function of N(0, 1). Then, each Zj has a marginal distribution Fj. The random vector Z is called a 

10 TGRF with symmetric positive definite (s.p.d.) dependence matrix St, denoted as TGRF n (.F, ^), 

11 where F — (Fx,..., F n ). The joint density of Z can be easily shown to be 

h(x) = (2tt)-5 |*|^exp (-^ T *^) n^T' ^ 

12 where /, is the density corresponding to Fj with parameters 0j, z = 1, . . . , n, is the density of 

13 N(0, 1), and e = JV^Fifo)}, . . . , ^{F^O}^. 

14 Clearly, a TGRF is obtained by transforming all the margins of a GRF with standard normal 

15 margins to desired marginal distributions Fj's. The resulting TGRF is not affected by the scales 
is of the original GRF since we can always standardize the margins to standard normals. The de- 
17 pendence structure of the TGRF, its copula, is still the Gaussian copula of the original GRF (e.g., 

Joe[ 1997, Nelsen, 2006). It is characterized by matrix but ^ no longer has the interpretation 



of correlation matrix. A TGRF n (.F, is completely specified by marginal distributions F and a 
Gaussian copula specified by a dispersion matrix More details and properties of the TGRF are 
presented in 201 1 University of Connecticut PhD thesis by M. O. Prates. 

22 A TGRF n (F, *) is a TGMRF if the GRF before the transformations is a standard GMRF with 

23 correlation matrix ^. As commonly used for GMRFs, it is more convenient to present a TGMRF 

24 using the precision matrix Q = \l> _1 , since it leads to an intuitive interpretation of conditional 

25 distributional properties. Let Z be a TGRF n (.F, Q^ 1 ) with a s.p.d. pre-transformation precision 

26 matrix Q. Since the transformations are marginal- wise, the Markov property is inherited by the 

27 TGMRF: for i j, _L Z^Zi^ if and only if = 0, where Zr^ is Z without the ith 

28 and jth observations. In other words, the precision matrix structure completely determines the 

29 conditional dependence structure of pairs given others. Since the undirected graph corresponding 

30 to a GMRF is retained in the resulting TGMRF, the equivalence of pairwise Markov property, local 



31 Markov property, and global Markov property for a GMRF (Ru e and Held[ |2005| ) are equivalent 

32 for a TGMRF. 

33 In the sequel, a TGMRF with marginal distributions F and precision matrix Q in the original 

34 GMRF scale is denoted as TGMRF n (i r ', Q). Matrix Q is not to be interpreted as precision but as 

35 a dependence matrix which characterizes the dependence structure. This property can be exploited 

36 in modeling practice to construct the precision matrix based on conditional dependences. 
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1 3 Spatial Generalized Linear Mixed Models 



2 The TGRF and TGMRF open a new avenue of random field models such as gamma field, beta field, 

3 and their Markov versions, which can be incorporated into the GLMM framework for modeling 

4 spatial dependence. Our departure point is the traditional GLMM with spatial random effects. 

5 Suppose that we observe (Yi, Xi) at sites i — 1, . . . , n, where Yi is the response variable and Xi 
e aqxl vector of covariates that correspond to response Yi at site i. Let e = (ei, . . . , e n ) T be a vector 
7 of unobserved random effects with joint distribution H, which introduces spatial dependence. A 
a spatial GLMM assumes that, given {X u e,), i = 1, . . . , n, the observations Yj's are independent 

9 with a distribution from the exponential family. Let /ij = E(Yi\X, e), where X = (Xi, . . . , X n ) T 

10 is the matrix of covariates. The conditional expectation /Xj is connected to the covariate X, t and 

11 random effect e« through a fixed link function g: 

g{fH) =Vi + (3) 

12 where r\i = Xj f3 is the fixed effect, and (3 is a q x 1 vector of regression coefficients of covariates 

13 Xi. The dependence among random effects e determines the spatial dependence among condi- 

14 tional means /j, = . . . , /i n ) T . Therefore, to fully specify a spatial GLMM, it is necessary to 

15 specify both the link function g and the joint distribution H of e. Commonly, H is chosen to be a 
is multivariate normal distribution with mean zero and covariance matrix S. 

17 Instead of introducing dependence among fi through the joint distribution H of random effects 

is e, we propose to specify a random field directly for fi. Specifically, our model for \x is 
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li ~ TGRF n (F, (4) 

where F = (Fi, . . . , F n ), Fj is the marginal distribution of fii, and ^ is the dispersion matrix 

20 characterizing the dependence structure of the underlying Gaussian copula. For independent data, 

21 in which case ^ is the identity matrix, this specification reduces to a class of GLMMs where 

22 the distribution of conditional mean /ij = + e^), instead of random effect e^, is specified. 

23 Our specification here is more general in that it incorporates dependence among all or part of /Vs 

24 through Gaussian copulas. 

25 The new model specifies the distribution of n through marginal distributions F and a Gaus- 

26 sian copula with dispersion matrix \&. It encompasses any model constructed from a link function 

27 g and H = N n (0, S) as a special case where F is the distribution function of /i, = g^ 1 ^ + e^), 

28 i — 1, . . . , n, and ^ is the correlation matrix of S. 

29 The TGRF model for fi provides a natural choice for the conditional means in hierarchical 

30 spatial models. For instances, one can use gamma margins for Poisson intensities and beta margins 

31 for Bernoulli rates, that can in turn be used, respectively, to model spatial count data or spatial 

32 binary data. The wide range of marginal distributions offer models that cover the traditional models 



33 as special cases and many more (e.g. |Prates et aL] |2010| ). The spatial dependence is completely 



34 characterized by the Gaussian copula, parameterized by the dispersion matrix ^ . For geostatistical 

35 modeling, where the observations sites may be irregularly spaced, one can parameterize ^ using, 



36 for examples, the exponential, spherical, or Matern structures (Banerjee et al. 2004 Ch.2). 
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Replacing the TGRF in model Q with a TGMRF, we model the conditional means fi by 

fi ~ TGMRF„(F, Q), (5) 

where the spatial dependence is characterized by Q, the precision matrix of the Gaussian copula. 
Since the copula is invariant to scale changes, we do not require that Q 1 is a correlation matrix 
as long as Q is scale free, s.p.d. precision matrix. 

Parameterization of Q is crucial and we propose to used the structure of the precision matrix 
of a conditional autoregressive (CAR) model ( |Besag[[T974 ). In a CAR model, the precision matrix 



is defined as Q/v, where Q determines the structure and v is a scale parameter. The scale v 
is not needed in our TGMRF model in ([5]). The structure Q is defined in such way that is 
nonzero if and only if site i and site j are neighbors of each other. To assure symmetry and positive 
definiteness, Q is defined as 

Q = M-\l-pW), (6) 

where M^ 1 is a diagonal matrix whose ith diagonal elements equal to the number of neigh- 
bors of site i, I is the identity matrix, p is a spatial dependence parameter, and W is a weight 
matrix providing contrasts of all neighbors to each site. Weight matrix W is determined by the 
neighboring structure and is of the form 



l/m, 

0, otherwise, 



15 where % ~ j indicates that site i is a neighbor of site j. 

is The proposed models fit naturally into the Bayesian framework. With carefully chosen priors 
17 for the parameters, Markov chain Monte Carlo (MCMC) algorithms can be developed to draw 



is samples from the posterior distribution of the parameters of interests (e.g., |Gelman et al.[ [2003 ). 

19 To compare different models for the same data, we propose to use the conditional predictive or- 

20 dinate (CPO) criterion (e.g., Gelfand et al. 1992[ Dey et al. 1997). The summary statistic is the 



21 logarithm of the pseudo-marginal likelihood (LPML), which is the summation of the log density of 

22 leave-one-out marginal posterior distribution. The performance of the CPO criterion in selecting 

23 the right models will be studied through simulations. The deviance information criterion (DIC) 

24 ([Spiegelhalter et al. 2002 ) is an alternative Bayesian model selection criterion. In our simulation 



25 studies, however, DIC had much higher variation than LPML and was outperformed in selecting 

26 the correct models. This might be explained by the fact that the DIC measures are highly depen- 

27 dent on the marginalization of the random effects, and become unstable when the distributions are 

28 nonnormal. 



29 4 Spatial Poisson Application 

30 Consider count data observed at n sites in a spatial domain. Let Yi be the count at site i, and with 

31 a q x 1 covariate vector JQ, i = 1, . . . , n. Poisson models are widely used for count data and the 

32 Poisson intensities are often modeled by gamma distributions. Few choices of gamma fields are 
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available in the literature. An exception is|Wolpert and Ickstadt|(|1998]), where a doubly stochastic 



2 process is used to construct positively autocorrelated intensity measures for spatial Poisson point 

3 processes which in turn are used to model the spatial count data. The TGMRF models provides 

4 new gamma Markov random fields to account for spatial dependence. 

5 4.1 TGMRF Models 

6 A GLMM introduces spatial dependence through a spatial random effect. Conditioning on = 

7 (/ii, . . . , /i„) T , the observed spatial count data Yj's are assumed to be independent, and each Yj is 
e Poisson with mean % = 1, . . . , n. The most commonly used GLMM for spatial count data uses 

9 the canonical log link on the Poisson intensities: 

log & = Xj/3 + ei, (7) 

10 where (3 is a q x 1 regression coefficient vector, e = (ei, . . . , e n ) T follows a GMRF with mean zero 

11 and a s.p.d. precision matrix ft/v, and v > is a parameter controlling the scale of the variance. 

12 Let of be the ith diagonal element of f2 _1 . Let F = (F\, . . . , F n ) T , where Fi is the distribution 

13 function of 

LN(X l T /3, vo\), z/ >0, i = l,..., n, (8) 

14 where LN(a, b) denotes a log-normal distribution with mean a and variance b on the log scale. It is 

15 clear that model ((T) is a special case of model (|5]) with Q = V l l 2 Q,V l l 2 and V = diag (erf, ... , of), 
is The TGMRF framework provides a new way to construct models for fi that incorporate spatial 
17 dependence and covariates. The Gaussian copula of TGMRFs captures the spatial dependence, 
is Any positive continuous distribution can be used to specify the marginal distribution of fi, and 

19 covariate effects can be accommodated into its parameters. Changing F in model (|5]) from log- 

20 normal to other distribution functions with positive support leads to new models. Gamma distri- 

21 bution is a natural choice for the margins. Let T(a, b) represent a gamma distribution with shape 

22 parameter a and scale parameter b, hence mean ab. Covariates can be incorporated into either one 

23 of the two parameters, resulting in two different gamma models as long as there is at least one co- 

24 variate. The gamma scale model, hereafter the GSC model, incorporates covariates into the scale 

25 parameter and defines the marginal distribution Fj as 

r(l/z/, v exp(Xj P)) , v>0, i = l,...,n. (9) 

26 The gamma shape model, hereafter the GSH model, incorporates covariates into the shape param- 

27 eter and defines the marginal distribution as 

T (exp(Xj/3)/u, u), u>0, i = l,...,n. (10) 

28 Under both models, the expectation of is the same, exp(Xj f3), but the parameter v has dif- 

29 ferent interpretations and should not be compared directly. TGMRF models with other marginal 

30 distribution for /ZjS can be constructed similarly. 

31 There is a subtle difference between the log-normal model ([8]), hereafter the LN model, and the 

32 two gamma models (|9]) and ( fTO] ). Unlike the gamma models, where the dependence structure does 
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1 not interfere with the marginal models, the dependence structure Q enters the marginal distribu- 

2 tions of /ijS through of in the LN model. This implies that a different model could be constructed 

3 with Fi being the distribution of 

LN(X7/3, u), u>0, i = l,...,n. (11) 

4 The variance parameter v could even incorporate covariates. These model could be used as alter- 

5 natives to the commonly used LN model ([8]) in the TGMRF framework. 



4.2 Abundance of Nenia tridens 
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Because of their abundance and critical roles in nutrient cycling, gastropods are of considerable 
ecological importance in terrestrial ecosystems (Mason} 1970[ ). In the Luquillo Mountains of 
Puerto Rico, Nenia tridens is one of the most abundant and widely distributed terrestrial gastropods 
in tabonuco forest flWfflig et al.[ [T998| |Bloch and WMg] |2006] |Willig et al.[ |20TTj ). Indeed, the 
forest ecosystems of the Luquillo Mountains have a long history of environmental study (e.g., 
Brown et al.[ 1983[ Reagan and Waide, 1996 ), resulting in deep understanding of the spatial and 



temporal dynamics of populations, communities, and biogeochemical processes, especially as they 
relate to natural and human disturbances ( Brokaw et al.[|20iTj ). 

Abundance data of N. tridens were collected from the Luquillo Forest Dynamics Plot (LFDP), 
a 16 hectare grid in tabonuco forest (18°20'N and 65°49'W), during the wet season of 1995 at 



each of 160 circular sites (3 m radius) on an lattice. As shown in Figure 1(a), there are 40 major 



sites in dark, 60 meters apart, and 120 supplementary sites in gray, 20 meters apart, placed inside 
the squares formed by the 40 major sites. Therefore, the data is available on a regular but sparse 
lattice. To define the graph for the TGMRF model, any two sites within 60 meters are considered 
neighbors, which results in different number of neighbors for major sites and for supplementary 



sites. Also shown in Figure 1(a) are an internal major site connected to its 20 neighbors and an 
internal supplementary site connected to its 16 neighbors. 

The abundance of N. tridens at each site was the minimum number known alive from four 



nocturnal surveys based on well established protocols on the LFDP ( |Willig et aL] 1998; Bloch and 
Willig., 2006). The observed count over the lattice is displayed in Figure [T(b) Possible covariates 



were topographic and habitat characteristics at each site. There were two topographic variables, 
elevation and slope. Four habitat variables were quantity of litter, canopy openness, apparency of 
sierra palm, and plant apparency. Quantity of litter was the mean number of leaves on the forest 
floor from each of four locations that were sampled at each site along mid-points of the radii from 
the center of the circle, arranged along cardinal compass directions (cardinal points). Canopy 
openness was the amount of light that penetrates to the understory (1.5 m above the forest floor) 
based on the mean number of open cross-hairs on a gridded densiometer, quantified from the four 
cardinal points. Plant apparency measured the volume of space in the understory that was occupied 
by plants using a plant apparency device at each of the four cardinal points, which captured the 
number of foliar intercepts along each of two perpendicular 1.0 m dowels placed at 0.5 m intervals 
from ground level to 3 m of height. Apparency of sierra palm measured specifically the apparency 
of Prestoea acuminata, a preferred substrate and food of N. tridens. 
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(a) 



(b) 



(c) 



Figure 1 : (a) Lattice of sampling sites at the LFDP in 1995 and the neighbor structure for an internal 
major site and an internal supplementary site, labeled by M and S, respectively, (b) Abundance of 
N. tridens. (c) Presence/absence of G. nigrolineata. 



1 We fitted Poisson regressions to the abundance data of N. tridens with three TGMRF models: 

2 LN, GSC, and GSH. For each model, the precision matrix Q of the Gaussian copula was specified 

3 with ([6]) from the CAR model. The prior distributions of regression coefficients % = 1, . . . , q, 

4 are independent N(0, 1/r), with r = 0.01 . The prior distribution of the scale or shape parameter v 

5 is specified as r(«i, k 2 ), with ki = 0.01 and n 2 = 100. These priors are set to be proper but vague 
b to allow the posterior estimates to be mainly data driven. Because we expected positive spatial 

7 dependence, a U (0, 1) prior is put on the spatial dependence parameter, p for the CAR model. 

8 The GSH model had the largest LPML, -482.12, followed by the GSC model (-483.90) and 

9 the LN model (—491.15). These results suggest that, The GSH model and the GSC model per- 

10 formed fairly closely, with the former being mildly preferred. The GSH model provides consid- 

11 erably better fit than the traditional LN model with a 9.1 difference in LPML. Since both models 

12 have the same number of parameters, the log Bayes factor is approximately twice LPML difference 
asymptotically (Gelfand and DeyJ 1994). From the rules suggested by Kass and Raftery ( ] 1995 ), a 
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14 log Bayes factor of 18.2, which falls in the category of 10 or higher, provides "very strong" evi- 

15 dence in favor of the GSH model over the LN model. As will be seen in our simulation study, when 
is the true model was the GSH model, 38 out of 100 replicates had LPML differences of greater than 
17 9.1 between the fitted GSH model and the fitted LN model; when the true model was the either one 
is of the two LN models considered, however, this rate became out of 100. 

19 The posterior point estimates and 95% highest posterior density (HPD) credible intervals of the 

20 parameters from the GSH model and the traditional LN model are summarized in Table [T] The two 

21 models lead to qualitatively the same conclusions. Neither elevation nor slope was found to have a 

22 significant effect on the abundance of N. tridens. Of the habitat variables, only canopy openness is 

23 negatively significant. More openness in the canopy implies fewer trees and dryer soil, which are 

24 not the preferred habitat condition by the N. tridens. The marginal scale parameter v is estimated 
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Table 1 : Posterior point estimates and 95% HPD credible intervals of the parameters in the Poisson 
regression for the abundance of N. tridens with the GSH model and the traditional LN model. The 
regression coefficients are in the order of intercept, elevation, slope, quantity of litter, canopy 
openness, plant apparency, and apparency of sierra palm. 

Parameters Specified Model 



GSH LN 





Estimate 


95% HPD 


Estimate 


f \ c frt t ttvt~\ 

95% HPD 


Regression coefficients 








A) 


1.990 


(0.841,2.731) 


2.137 


(1.434,2.839) 


fa 


-0.062 


(-0.293,0.174) 


-0.037 


(-0.329, 0.226) 


02 


-0.046 


(-0.179, 0.081) 


-0.045 


(-0.163,0.064) 


03 


-0.103 


(-0.238, 0.020) 


-0.105 


(-0.232, 0.015) 


04 


-0.143 


(-0.292, -0.014) 


-0.129 


(-0.252, 0.007) 


05 


0.027 


(-0.100, 0.161) 


0.021 


(-0.107,0.139) 


06 


0.033 


(-0.104, 0.161) 


0.024 


(-0.105,0.149) 


Scale and spatial dependence parameters 






V 


4.924 


(3.902,7.312) 


5.134 


(3.534, 6.861) 


X 


0.951 


(0.851,0.996) 


0.953 


(0.851,0.999) 



to be 4.924. The spatial dependence parameter p is estimated as 0.951, with a HPD interval away 
from zero, which indicates a higher spatial dependence in the model is needed. 

4.3 Simulation Study 

To assess the fitting capacity of the TGMRF models, the properties of the Bayesian inferences, 
and the effectiveness of LPML as a model comparison criterion in this context, we conducted a 



b simulation study using the lattice and neighbor structure in Figure 1(a) Each of the three models 

7 was used as data generating models. In addition to the intercept, one covariate was generated from 

8 iV(0, 1), and the true covariate coefficient vector was = (1.0, 0.7). The precision matrix of the 

9 TGMRF took the form of (|6]) for the CAR model, with p = 0.8. The parameter v, which is related 

10 to the variance in all models, was set at v — 2, although it has completely different meanings. 

11 With v = 2, the gamma scale model and the gamma shape model appeared to be more similar 

12 to each other than to the log-normal model. To make a more interesting comparison, a second 

13 log-normal model was also used to generate data, where v = 6.5 was chosen because it provides 

14 good approximation to the gamma scale model with v = 2. In summary, we had a total of four 

15 data generating models: two LN models LN1 and LN2, one GSC model, and one GSH model. 

is For each data generating model, we generated 100 datasets, and fit each dataset with all three 

17 proposed TGMRF models. In each fitting process, a vague prior, r(0.01, 100), was set for the 

is dispersion parameter v, and an uninformative U (0, 1) prior was set for the spatial dependence pa- 

19 rameter p. Independent iV(0, 100) priors were set on regression coefficients 0. Table [2] summarizes 
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Table 2: Summaries of posterior mean, standard deviations (SD), and LPML from 100 replicates 
in the simulation of spatial Poisson regression. 



True 
IVTdHpI 


Param 


True 

\/a 1 1 l p 
V ell u L. 






Specified Model 






T N 

J— /I > 








GSH 

VJ Oil 




IVlCdll 


en 
o u 


Mean 


SD 


ivicdn 


o u 


LN1 


A) 


1.00 


0.99 


0.10 


1.08 


0.09 


l.n 


0.14 




fa 


0.70 


0.70 


0.06 


0.70 


0.05 


0.67 


0.06 




p 


0.80 


0.53 


0.26 


0.55 


0.26 


0.55 


0.26 




V 


2.00 


2.32 


0.73 


6.35 


1.53 


0.84 


0.46 




LPML 




-331.90 


10.76 


-332.55 


10.87 


-335.57 


11.11 


LN2 


A) 


1.00 


0.98 


0.16 


1.33 


0.18 


1.46 


0.23 




fa 


0.70 


0.70 


0.08 


0.70 


0.07 


0.58 


0.08 




p 


0.80 


0.60 


0.22 


0.61 


0.21 


0.64 


0.22 




V 


6.50 


6.97 


1.42 


2.00 


0.41 


3.50 


1.16 




LPML 




-362.90 


15.06 


-365.69 


14.78 


-371.35 


15.77 


GSC 


0o 


1.00 


0.76 


0.18 


0.99 


0.13 


1.09 


0.19 




0i 


0.70 


0.70 


0.08 


0.70 


0.07 


0.61 


0.07 




P 


0.80 


0.59 


0.23 


0.59 


0.23 


0.60 


0.23 




V 


2.00 


6.34 


1.30 


2.24 


0.53 


2.18 


0.73 




LPML 




-336.34 


14.96 


-335.38 


14.47 


-340.95 


15.12 


GSH 


0o 


1.00 


0.71 


0.20 


1.00 


0.14 


0.99 


0.18 




0i 


0.70 


0.77 


0.08 


0.71 


0.08 


0.70 


0.07 




P 


0.80 


0.64 


0.21 


0.62 


0.21 


0.63 


0.21 




V 


2.00 


7.09 


1.31 


1.95 


0.48 


2.17 


0.60 




LPML 




-335.86 


17.27 


-335.96 


16.60 


-328.81 


17.35 



the mean and standard deviations of the Bayesian estimate of the parameters and LPML from the 
100 replicates. 

When the model was correctly specified, the true values of the regression coefficients were 
recovered very well. The estimates seems to be upward biased for the dispersion parameter v but 
downward biased for the dependence parameter p, suggesting that spatial dependence and spatial 
heterogeneity are hard to identify. When the model was misspecified, the regression coefficient 
estimates were still recovered reasonably well, especially in the GSC model and the GSH model, 
probably because the mean of fi was still correctly specified, regardless of the misspecified model. 
In all cases, the average of the LPML statistic was higher for correctly specified models than for 
the misspecified models, with similar variation under different models. 

To gain a clearer picture on model comparison using LPML, we summarize the frequencies of 
the models selected with the highest LPML from all 100 replicates under each of the four models 
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Table 3: Frequencies of selected model using the LPML statistics for the 100 simulated datasets. 

True model Frequency selected 





LN 


GSC 


GSH 


LN1 


59 


29 


12 


LN2 


77 


16 


7 


GSC 


34 


59 


7 


GSH 


6 


5 


89 
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True model 

(a) 



GSH 



log it 1 



logit 2 
True model 

(b) 



beta-log it 



Figure 
ulation 



2: LPML difference between the correct model and misspecified models, (a) Poisson sim- 
. (b) Bernoulli simulation. 



(Table [3]). The criterion seems to be very effective when the true model was the GSH model, 
correctly selecting the true model 89 times. When the true model was LN1 or GSC, the correct 
model was selected 59 times in either case with our sample size, while the alternative GSC model 
or LN model was selected 29 and 34 times, respectively; the GSH model was selected only 12 
and 7 times, respectively. This indicates that the LN model and the GSC model provides good 
approximation to each other, similar to their well known similarity in univariate modeling without 
covariates and spatial concerns; a large sample would be necessary to distinguish them effectively. 
With our sample size, when the true model was LN2, the LPML was able to differentiate the 
LN model better from the GSC model, correctly selecting the LN model 77 times. Therefore, 
the similarity between the GSC model and the LN model appear to be different under different 
scenarios. The GSH model seems to have specific characteristics that make it further away from 
both the LN model and the GSC model in the model space. 
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A closer look at the difference in LPML across models is through box plots. Figure |2(a) 



presents the box plots of the difference in LPML between the correct model and two misspecified 
models for each true model. The magnitude of the differences provides guidance in practice on 
what models are similar to each other and on how big a difference is important. In the spatial 
setting we considered, the LN model and the GSC model were very similar, as seen from the boxes 
centered near zero. The majority of each box plot is well above —5, suggesting that if the LPML 
of one model is observed to be higher than that of another model by 5, then it is very unlikely that 
the other model is the true model. 



5 Spatial Bernoulli Application 

Consider presence/absence data at n sites in a spatial domain. Let Yi be 1 if presence is observed 
and otherwise at site i, with a q x 1 covariate vector X u i — 1, . . . , n. 



2 5.1 TGMRF Models 

13 Conditioning on fi = (hi, . . . ,/i n ) T , the observed data Fj's are assumed to be independent, and 
u each Yi is Bernoulli with mean fa, % — 1, . . . , n. The traditional spatial GLMM for binary data is 

T, 



logitOO = Xl /3 + e. 
15 where (3 is a q x 1 regression coefficient vector, e = (ei, 



(12) 

e n ) T follows a GMRF with mean 

is zero and precision matrix ft/u, and v > is a parameter controlling the scale of the variance. 

17 Let of be the ith diagonal element of Q _1 . Let F = . . . , F n ) T , where F { is the distribution 

18 function /i« = logit -1 (X 4 T /3 + i = 1, . . . , n. Then, model ( |T2| ) is a special case of model ([5]) 

19 with Q = V 1/2 ttV 1/2 and V = diag (of, . . . , a*). 

20 Changing F in model ([5]) to any distribution function defined over the (0, 1) support leads 

21 to new models. Covariate effects can be accommodated into the marginal parameters. Spatial 

22 dependence is modeled through the Gaussian copula with dispersion matrix fl^ 1 . 

23 The beta distribution is a natural choice for the margins. Let Beta(z/p, v(l — p)) represent a beta 

24 distribution with mean parameter p and dispersion parameter v. Covariates can be incorporated into 

25 Jthe mean parameter p using any transformation function from 3? to (0,1) (e.g., Ferrari and CribarT] 

26 Neto[ |2004|). We propose a beta-logit model that incorporates covariates into the mean parameter 



27 p using a inverse logit transformation and defines marginal distribution F± as 



Beta 



exp(XTff) 



exp(X t T /3) 



U exp(Xj(3) + 1 ' V \ exp(XjP) + 1 



v > 0, 



n. 



(13) 



There is again a subtle difference between the logit model ([12]) in comparison with the beta 



28 

29 logit model ( fT3[ ). In the logit model ( fT2| ), parameters in the dependence structure Q enters the 

30 marginal distribution, whereas it does not do so in the beta model. 
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5.2 Presence of Gaeotis nigrolineata 
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19 



25 



26 



28 



29 



30 



36 



37 



38 



39 



40 



Gaeotis nigrolineata is a common terrestrial gastropod in tabonuco forest of the Luquillo Moun- 
tains of Puerto Rico ( |Willig et al. [ [T998| |Bloch and WmrgTj [20061 |Willig et al.[ |20TT] ). Its spatial 
distributions is believed to be associated with the abundance of live sierra palms, its preferred sub- 
strate. Generally, Gaeotis nigrolineata is less abundant than is N. tridens. It often occurs in low 
numbers, and is characteristically absent from a significant proportion of the sites across the LFDP. 
Therefore, it is more suitable to analyze the presence/absence data for this taxon. 

The presence/absence data were obtained by dichotomizing the abundance of G. nigrolineata, 
which were determined in the same manner as described for N. tridens (Section [472| ). In particular, 
we have one for presence and zero for absence at each site. The distribution of incidences for G. 
nigrolineata is apparently heterogeneous with spatial clustering across the data collection lattice; 
see Figure 1(c) All but one of the covariates as described in Section 4.2 were used to model spatial 
dynamics of G. nigrolineata. Since G. nigrolineata does not live or feed in the leaf litter, quantity 
of litter was not included as a covariate in its analysis. 

We fitted Bernoulli regressions for presence/absence data of G. nigrolineata with two TGMRF 
models: logit and beta-logit with precision matrix of the CAR model. Prior distributions for the 



models parameters were selected the same as those described in Section 4.2 

The LPML values were —99.31 and —103.51 for the beta-logit model and the logit model, 
respectively. Therefore, using a CAR dependence structure, the beta-logit model fits better than 
the traditional logit model. The approximate log Bayes factor was 8.4, which falls in the category 
of [6, 10) suggested by |Kass and Raftery (1995]), "strong" evidence favoring the beta-logit model 
over the logit model. To be seen in our simulation study, when the true model was beta-logit, 68 
out of 100 replicates had LPML differences of greater than 4.2 between the fitted beta-logit model 
and the fitted logit model, but the rate was 1 or out 100 when the true model was a logit model. 

The posterior point estimates and 95% HPD credible intervals for parameters in both models 
are summarized in Table |4j The conclusions of the two models are virtually the same. Neither 
elevation nor slope had a significant effect on the incidence of G. nigrolineata, as in the case for N. 
tridens. Of the habitat characteristics, only plant apparency had a significantly negative effect on 
the incidence of G. nigrolineata, That is, the greater the volume of vegetation in the understory of 
the forest, the lower the abundance of G. nigrolineata. The apparency of sierra palm, which mea- 
sures the preferred substrate for the G. nigrolineata, was found to be almost positively significant 
with the 95% HPD credible interval barely including zero. The negative effect of plant apparency 
was surprising but the paradox may be resolved if high plant apparency in the understory indicates 
the presence of an opening in the canopy, and attendant temperatures (high) and humidities (low) 
outside of the fundamental niche of G. nigrolineata, precluding its presence even though its pre- 
ferred substrate may be common. The spatial dependence parameter p is estimated as 0.760 and 
0.803 in the two models, respectively, indicating strong spatial dependence within neighbors areas. 

It is worth noting that although the beta-logit model agrees with the logit model in the direc- 
tions of the covariates effects, it has much smaller widths in the HPD credible interval does the 
logit model. This indicates better precision of the estimating the coefficients. The v parameter 
does not have the same interpretation in the two models, and, hence, they are not directly compa- 
rable. Nevertheless, from Table |4} we can see that for the beta-logit model where v is a marginal 
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Table 4: Posterior point estimates and 95% HPD credible intervals of the parameters in the 
Bernoulli regressions for the presence/absence of G. nigrolineata with the beta-logit model and 
the traditional logit model using the CAR dependence structure. The regression coefficients are in 
the order of intercept, elevation, slope, canopy openness, plant apparency, and apparency of sierra 
palm. 

Parameters Specified Model 

beta logit logit 
Estimate 95% HPD Estimate 95% HPD 



Regression coefficients 



Pa 


0.298 


(-0.481, 1.253) 


0.226 


(-1.733, 1.839) 


Pi 


0.326 


(-0.174, 0.777) 


0.527 


(-0.578, 1.744) 


fa 


0.087 


(-0.238, 0.428) 


0.134 


(-0.521,0.803) 


Ps 


-0.014 


(-0.340, 0.344) 


-0.051 


(-0.715,0.704) 


Pa 


-0.500 


(-0.887, -0.137) 


-0.894 


(-1.784, -0.174) 


Pb 


0.270 


(-0.074, 0.644) 


0.538 


(-0.200, 1.330) 


e and spatial dependence parameter 






V 


1.699 


(0.235, 4.551) 


92.808 


(0.200, 234.300) 


P 


0.760 


(0.264, 0.998) 


0.803 


(0.326, 0.999) 



overdispersion parameter, v is more identifiable with a small HPD interval. For the logit model, 
the v parameter is the marginal variance. The estimate implies a standard deviation of 9.634 with a 
wide 95% credible interval of (0.447, 15.297). On the log scale of fj,, such a magnitude of variation 
may not mean much on the original scale of fj, since the log transformation explodes at zero, and 
this may explain the poor identification of the spatial logit model. 



5.3 Simulation Study 

A simulation study was conducted for the spatial Bernoulli regressions. Both the logit model and 
the beta-logit model with the CAR dependence structure were used to generate data. Except for the 



response variable, the simulation setup was the same as that in Section 4.3 with model parameters 
P = (1.0, 0.7), p = 0.8 and v — 2. Again, since v has different interpretation in the two models, 
a second logit model with v — 1 was also used to generate data in attempt to approximate the 
beta-logit model with v — 2. For each of three true models, we generated 100 datasets, and fit 
each dataset with each of two TGMRF models. The priors were chosen in the same manner as 



Section 4.3 Table |5] summarizes the posterior mean and standard deviations estimates from 100 



replicates. 

Similar to the results from Section 4.3 when the model was specified correctly, the true values 
of regression coefficients are recovered very well; the dispersion parameter estimate tended to be 
bigger than true value; and the dependence parameter estimate appeared to be downward biased. 
When the true model was the beta-logit model, the average LPML value of the beta-logit model 
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Table 5: Summaries of posterior mean, standard deviations (SD) and LPML from 100 replicates in 
the simulation of spatial Bernoulli regression. 



True 


Par am 


True 


Specified Model 




Model 




Value 


lopit 




beta-logit 








Mean 


SD 


Mean 


SD 


logit 1 


A> 


1.00 


1.02 


0.22 


0.95 


0.23 






0.70 


0.72 


0.23 


0.65 


0.20 




p 


0.80 


0.50 


0.29 


0.47 


0.27 




V 


2.00 


2.33 


6.03 


4.24 


2.49 




LPML 




-92.30 


5.59 


-91.77 


5.77 


logit 2 


A> 


1.00 


1.03 


0.22 


0.97 


0.23 




fa 


0.70 


0.73 


0.22 


0.66 


0.20 




p 


0.80 


0.50 


0.29 


0.46 


0.27 




V 


1.00 


1.56 


3.57 


3.79 


2.52 




LPML 




-91.67 


5.53 


-91.31 


5.58 


beta-logit 


0o 


1.08 


0.22 


0.18 


1.01 


0.27 




0i 


0.70 


0.78 


0.23 


0.68 


0.20 




P 


0.80 


0.52 


0.29 


0.56 


0.27 




V 


2.00 


1.10 


1.09 


3.99 


2.37 




LPML 




-96.16 


6.27 


-88.16 


8.21 
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Table 6: Frequencies of model selection using the LPML statistics for the 100 simulated datasets. 



True model 


Frequency Selected 


logit beta-logit 


logit 1 


46 54 


logit 2 


49 51 


beta-logit 


16 84 



1 was 8 higher than that of the logit model. When the true model was the logit 1 or logit 2, however, 

2 the average LPML value of the beta-logit model was very close to (actually slightly higher than) 

3 that of the logit model in both cases. This implies that the beta-logit model is quite accommodating 

4 and can provide close approximation to the logit model; with the sample size in our simulation, 

5 they are hard to distinguish. 

6 Table [6] summarizes the frequencies of the models selected with the highest LPML from all 

7 100 datasets generated under each scenario. When the true model was the beta-logit model, the 
s LPML criterion worked effectively, correctly selecting the true model 84 times. When the true 
9 model was logit 1 or logit 2, however, the logit model and the beta-logit model were selected with 

almost equal frequency, indicating that the beta-logit model provides very good approximation of 

11 the logit model with our sample size. 

12 Box plots of the difference in LPML between the correct model and the misspecified model 
are shown in Figure |2(b) The boxes are surprisingly tight around zero when the true model is 
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13 



14 the logit model, indicating that the beta-logit model approximates the logit model very closely in 

15 terms of LPML. When the true model was the beta-logit model, however, the LPML value of the 
is logit model was very unlikely to be higher than that of the correctly specified model. The majority 
17 of all box plots were well above —5. A difference of 4.2 between the two models as observed in 
is the analysis of presence/absence of G. nigrolineata seems to be quite strong evidence in favor of 
19 the beta-logit model. 



20 



6 Discussion 



23 



30 



In geo statistics, the trans-Gaussian kriging approach if often used to transform the responses to 
achieve joint normality ( [Cressie |1993| ). Although the dependence structure in a trans-Gaussian 
kriging approach is also a Gaussian copula, our approach is different in several aspects. Our trans- 
formation is not to Gaussian but from Gaussian, and our model is directly built for the variable 
of interest, rather than on some power transformation of it, which may be hard to interpret. Even 
when viewed as a to-Gaussian transformation, our transformation is margin specific and can incor- 
porate covariates. From a hierarchical model point of view, our random fields are mostly useful 
for model parameters such as Poisson intensity or Bernoulli rate, which is a different domain than 



29 kriging or spatial interpolation, (see, De Oliveira et al. 1997; Azzalini and Capitanio 1999) 



The proposed models are highly likely to be favored by the LPML model selection criterion 
when they are the true models. Even when they are misspecified, they may still be competitive 
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1 by providing a close approximation in the misspecified class for small to moderate sample sizes 

2 in practice. Our simulation study for the spatial Poisson regression examined the performances 

3 of three TGMRF models with different marginal distributions and different parameterizations to 

4 incorporate covariates. The GSC model appears to be more versatile than the traditional LN model 

5 in that, even when the latter is the true model, the former may provide a very close approximation 
e under practical sample size. The gamma shape model provides another way to improve data fitting. 
7 For the abundance date of N. tridens, the GSH model provided the best fit and sheds light on 
s important predictors. In the simulation of the spatial Bernoulli regression, the beta-logit model 

9 appeared to be as good as the logit model in terms LPML even when the data were generated by 

10 the logit model, but the opposite was not true. For the presence of G. nigrolineata, the beta-logit 

11 improved the fitting with narrower HPD credible intervals. In real world applications, where the 

12 true model is unknown, the class of our proposed models may be useful in approximating the 

13 unknown truth. 
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